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We study the Navier-Stokes steady states for a low density monodisperse hard sphere 
granular gas (i.e, a hard sphere ideal monatomic gas with inelastic interparticle colli- 
sions) . We present a classification of the uniform steady states that can arise from shear 
and temperature (or energy input) applied at the boundaries (parallel walls). We con- 
sider both symmetric and asymmetric boundary conditions and find steady states not 
previously reported, including sheared states with linear temperature profiles. We pro- 
vide explicit expressions for the hydrodynamic profiles for all these steady states. Our 
results are validated by the numerical solution of the Boltzmann kinetic equation for the 
granular gas obtained by the direct simulation Monte Carlo method, and by molecular 
dynamics simulations. We discuss the physical origin of the new steady states and derive 
conditions for the validity of Navier-Stokes hydrodynamics. 



1. Introduction 

The kinetic theory of non-uniform gases has advanced significantly in the last few 
decades ( |Hilbert|1912| [Chapman fc Cowling|1970||Grad|1963||Ernst|1973[ ). For instance, 
the Chapman-Enskog method (Chapman & Cowling 19701, which allows the determi- 



nation of the transport coefficients associated with the equations of evolution of the 
average fields n, v and T (density, flow velocity and granular temperature, respectively), 
has recently been successfully extended to granular gases, where energy is not conserved 
in the collisions (Goldhirsch [2003). However, debate continues regarding fundamental 



issues such as the validity, in granular gases, of the assumptions made in the standard 
Chapman-Enskog method ( Goldhirsch[[2003 ) . Nonetheless, this work has produced a se- 
ries of results that can help explain important features observed in real or computer 



experiments (for example, non-equipartition in granular gases mixtures, Garzo & Dufty 



1999a[[Barrat fc Trizac[[2002a[ [Feitosa fc Menon[[2002[ [Wildman fc Parker[[2002 1 

The use, in the macroscopic balance equations for a granular gas, of the constitutive 
equations of the stress tensor and heat flux and the corresponding hydrodynamic coeffi- 
cients yields a peculiar hydrodynamics from which many interesting transport phenomena 
arise ( Goldhirsch[[2003 ) . In the hydrodynamic equations a nonuniform temperature pro- 
file arises in general from the effects of viscous heating (in cases where there is shear), 
coUisional cooling, and thermal transport. It has been shown that a hard-disk or hard- 
sphere system that is sheared and heated by two identical infinite parallel walls can have a 
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temperature profile of negative curvature (when viscous heating exceeds collisional cool- 
ing, s ee works byjJenkins fc Savage„1983] |Nott, Alam, Agarwal, Jackson fc Sundaresan| 



1999 Tij, Tahiri, Montanero, Garzo, Santos fc Dufty|2001 l, positive curvature (when the 
reverse is true), or a flat temperature profile ( Campbell|1989 |. This last condition, occur- 
ring when viscous heating and collisional cooling exactly balance, is called Uniform Shear 
Flow (USF), and has been extensively studied ( Goldhirsch|2003 1 . In the case of a heated 
gas (with no shear), temperature profiles with positive curvature ( Barrat fc Trizac|2002& I 
have been reported. These studies considered only symmetric boundary conditions for 
the granular temperature (for example, both walls at the same temperature). Alam, Nott 
and collaborators (see, for instance 'Alam fc Nott|1998 Nott et al.\1999\ and Wang et al. 



{ 1996| ) have obtained stability and bifurcation criteria for the base steady flows resulting 
from symmetric boundary conditions for the temperature, with specific studies on the 
stability of the uniform shear fiow ( Ala m|2006 1 . On the other hand, there are also works 
where asymmetric boundary conditions are analyzed. In these cases, except for the work 



by Brey et al. (2001), where a linear temperature profile in an unsheared granular gas 



is found, the same types of profiles as in the symmetric case are reported: for instance. 



Khain et al. ( 2008 ) derive theoretical criteria for the correction of the temperature jump 



at the boundary, and Galvin et al. ( 2007 ) study via MD simulations the deviations from 



NS behaviour in the boundary layer. 

In this work we focus on the systematic derivation of all steady profiles and we show 
that the number of steady states arising from asymmetric boundary conditions for the 
granular temperature is actually twice the number of types of profiles reported previously. 
With 'type of profile' we mean profiles that share basic properties like a given temperature 
curvature, a particular form of the velocity profile and a specific balance between viscous 
heating, collisional cooling, and the exchange of heat fiow with the boundaries. We also 
describe in this work the basic balance mechanisms of heat flow, viscous heating and 
collisional cooling that can occur in the steady base states of a granular gas. In effect, 
we demonstrate that the temperature profile curvature (positive, negative or zero) is 
determined by an interplay of collisional cooling and viscous heating in the bulk of the 
fluid and also the difference between the granular temperature at the different boundaries 
of the system. This produces some surprising results, for instance, we find that effects 
of the temperature boundary conditions may compensate collisional cooling in the bulk 
even in the absence of shearing and produce a steady temperature profile with negative 
curvature and we find that it is possible for a sheared granular gas to have a perfectly 
linear temperature profile and that this happens even though collisional cooling and 
viscous heating are not balanced (contrary to what happens in the USF). We validate 
our resulats with simulation data (both Monte Carlo method and molecular dynamics). 

The structure of the paper is as follows: In §|2]we determine the differential equations 
for the steady temperature and flow velocity proflles that result from the Navier-Stokes 
equations for the granular gas. In § |3l we obtain all possible types of unidimensional 
steady profiles (henceforth, steady base states): in § 3.1 we write the differential equations 

n 

we 



in reduced magnitudes; the resulting form of these equations allows us to define in 
the relevant reduced (microscopic/hydrodynamic) length and time scales; in 



3.3 



obtain analytical solutions of the steady hydrodynamic profiles and describe the general 
properties of the resulting temperature profiles; finally, in § |3.4| we give the classification 
of the different types of base steady states and flows in granular gases, in the case of 
constant pressure. We end with § El where we compare our theory with simulation data 
obtained by the direct simulation Monte Carlo method (DSMC) and molecular dynamics 
simulations (MD). 
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2. Continuum balance equations and transport coefficients 

With the assumptions of mass and momentum conservation in the coUisions and from 
integration of the Boltzmann equation for a granular gas x{l,mv,mv^/2} the mass, 



momentum and energy balance equations become ( Jenkins fc Richman||1985| Brey et al. 
1998] |Sela fc Goldhirsch||1998D 



Dn 



-nV • u 



(2.1) 



Du 

lot 



V-P 



mn 



(2.2) 



DT 
~Dt 



TC = -—{P 
an 



V u+V -q) 



(2.3) 



where D/Dt = d/dt + it • V is the material derivative, m the particle mass, and n is 
the particle density, u the flow velocity, T the temperature, P the stress tensor and q 
the heat flux (see [Brey et al.\1998 Brey fc Cubero|2001 Goldhirsch|2003 and references 
therein for more details on the deflnitions and derivations). The rate at which kinetic 
energy is lost in the collisions, the cooling rate, is denoted as (. The balance in (2.3 1 



between the term coming from viscous heating (the first term on the right hand side) 
and the coUisional cooling (the second term on the left hand side) will determine the 
gradient of the heat flux in the steady state {d/dt = 0). 

The hydrodynamic description (and thus, our results) is valid only for sufficiently small 
gradients of n, u, T, meaning that the typical microscopic length scale A (the particle 
mean free path) must be small compared to a typical system macroscopic length scale L 
(determined by the gradients of n, u, T): X/L ^ 1. 

For small enough X/L, we may retain in the Chapman-Enskog expansion only terms 
up to first order in the gradients ( Chapman fc Cowling|[l970 ) , and in this case we find 
the following constitutive relations for the stress tensor and heat flux ( Brey et ar]|1998 ) 



P — pi — rj 



Vu 



Vm^--(V-m)/ 



(2.4) 



q= -kVT- ^iVn, (2.5) 

where p = nT is the hydrostatic pressure, the coefficients rj and k are known as the 
viscosity and heat conductivity respectively, and /i is a coefficient that only appears when 
energy is not conserved in the collisions ( Brey et aZ.||1998 ). Inserting (2.4|, (2.51 into the 
balance equations (2.1 )-(2.3) we obtain the Navier-Stokes equations for the granular gas. 



The transport coefficients for a low density granular gas of hard spheres have been cal- 



culated by Sela fc Goldhirsch (19981 and Brey et al. (1998). We use here the expressions 



for hard spheres (with the same notation) calculated by Brey et al. ( 1998 1 (see the more 
recent work, jBrey fc Cubero|[2001[ for more general expressions for a d-dimensional sys- 
tem). From dimensional considerations in ( |2.1| -(2.3l with (2.4l-(2.5l, the hydrodynamic 
Navier-Stokes coefficients have the following form (Brey et al. 19981 



1/2 



where tjq 






.SrV2 



7-3/2 



M = Mo- 



C = Co; 



(2.6) 



rpl/2 ' 

= KQK*{a), /ij5 = Ko/i*(a) and Q = C*i'^)/Vo ^^^ constants, 
are related respectively to the viscosity, r]o{T) = tj'qT^^'^, and 
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Figure 1. Schematic view of the system subject of study. The granular gas is heated and 
sheared from two infinite parallel walls. 



thermal conductivity, Kq{T) = kIi^T^/"^, of an elastic hard disk/sphere gas (Chapman &; 
[Cowling,,1970i |Brey e^ ar]|1998| ), with ( [Chapman fc CowUng||l97n| : 

V', = ^r(d/2)7r 



- 1/2 - 



(d-i) 



d{d+2f 
16(d-i; 



r(d/2)7r" 



t^-l/2^-(d-l) 



(2.7) 



The dimensionless functions 7f{a), fi*(a), K*(a) and C*('^) depend on the coefficient 
of normal restitution a and in the elastic limit ?7*(1), «*(!) = 1 and /^*(1), C*(l) = 0- 

Transport coefficients in th e quasielastic limit (|Sela fc Goldhirsch||1998|) and for mod- 
erate density granular gases (Jenkins fc Richman 1985| [Garzo fc Dufty 1999 &) and co- 
efficients obtained from expansions from arbitrary reference states ( Lutsko||2006 Garzo 
12006) can also be used. Also, alternate hydrodynamic descriptions based on BGK-type 



kinetic models (Vega Reyes, Garzo fc Santos 20071 and coefficients calculated with a 
more refined Sonine expansions (Garzo, Santos fc Montanero 2007 Garzo, Vega Reyes 
fc Montanero||2009 1 may be used. 



2.1. Hydrodynamic steady states 

Let us consider first a system with two parallel walls in the planes y = —h/2 and y = 
+h/2. The system is infinite in the x and z directions. The walls are kept at temperatures 
Tyj{y = +h/2) = T^u and Ti,{y = -h/2) = T^d- Additionally, the walls at y = -h/2 
and y = +h/2 may be moving in the x direction, in general with different velocities, 
so that the granular fluid adjacent to the walls has velocities v^{—h/2) — Uum^x and 
Vw{+h/2) = Uwu^x respectively, where Bj is a unit vector in the j direction. A schematic 
view of the system is depicted in figure [l] We will study steady states {d/dt = 0) that 
are uniform in the direction of the flow (dux/dx = 0), and assume that there is nothing 
to prevent dp/dx = (a gravitational force in that direction, for instance) and we 
will also assume that the system is uniform in the z direction, d/dz = 0. (Thus all 
conclusions below also apply to a system of sheared two-dimensional inelastic disks, 
removing coordinate z). Under these conditions, our steady states will have the general 
form: u — u{y)ex, T — T(y), while p is a constant. 

From ( |2.2[ ) and ( |2.3[ ), and taking into account the system symmetries, we derive the 
Navier-Stokes equations for the steady states, 



Trlk) 



1_ 



(2.8) 
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d C P^ 






dux 
dy 



(2.9) 



Kq/3*(q!) and we multi plied niomentum and energy 



T/Tr respectively. In (2.8 1 we have integrated once 



where /3o — '«o('**(q^) ~ M*(q^)) = 
balance equations by l/-\/7V and 
over y, with 7 the resulting integration constant, which we call shear rate; and in (2.9 1 
we have taken into account that for constant pressure {T^^'^ /n){dn/dy) — —VT{dT/dy). 
Thus, the hydrodynamic problem needs 4 boundary conditions to be closed (two for 
each T and Ux)- We will not specify a relation between the hydrodynamic fields at 
the boundaries and the boundary conditions themselves (i.e., any functions u{vtu,Tu,), 
T{Viu,T^) at the boundaries). In ^3.3 we discuss the kinds of physical boundaries that 
would produce the specific profiles we describe. 



3. Analytic solution of unidimensional steady states Avith zero 
pressure gradient 

3.1. Reduced differential equations 

For the microscopic length and time scales we will use the mean free path A and collision 
frequency v defined at an arbitrary reference point, A^ and Vr, determined according to 



Xr ^ AdiV2nra^'^-^^y 



(3.1) 



2Tr 



^d-1 



A, 



(3.2) 



In (3.1 ), Ajj — ^/2'^^T{d/2)TT 2- is a constant that depends on the system's dimen- 
sion d, and Tr and n^ = Pr/Tr indicate respectively the values of temperature and density 
at an (initially arbitrary) reference poin t. 



All of the partial derivatives in (2.8 1 and (2.9) have a prefactor \/T, so that with 
the use of a scaled variable / such that y^T/Trd/dy — d/dl, the previously non- linear 
differential equations (2.8) and (2.9 1 become linear in the new variable (a similar change 



of variable was first used by [Santos, Brey fc Garzo 1986 for studying the steady heat 
flow in an gas of elastic spheres). The use of the extra factor \J\jTr allows us to preserve 
the dimensionality of the spatial gradients in the new variable. Notice also that the 
resulting summands in (2.8) and (2.9) have dimensions of v^ and vnv\ respectively; i.e., 
they set two characteristic (inverse) hydrodynamic time scales for the gradients of it 
and T. Thus, expressing equations (2.8), (2.9) as a function of derivatives in the scaled 
variable I and dividing them by v^ and mv^ respectively, we obtain the following linear 
and dimensionless differential equations 



dUx 

di 



rj*{a) 



Via) 



(3.3) 



(3.4) 



where we denote the corresponding dimensionless magnitude of T as T = T/T^, y as 
y = y/Xr, Ux as Ux = Ux/{XrVr) etc. 
The parameters 7 and T{a) read 
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7 
P 



(3.5) 



r(a) = (j)d 



1 



(3*{a)r]*{a) 



-C{aW{a)~f 



(3.6) 



d-l 



where <^d--dp+2) ■ 

The dimensionless magnitude r(Q;) represents the reduced coUisional cooling rate minus 
the viscous heating rate (coming from the term proportional to 7^), and because of this 
we will call it the effective cooling rate. The equilibrium between collisional cooling and 
viscous heating determines the curvature of the temperature as a function of the scaled 
variable I, but, as we will see in § 3.3 it does not the determine the curvature of Tijj). 



3.2. Reduced time and length scales 



The solution of (3.4 1, valid for all values of 7 and r(a), in terms of the variable /, is 
straightforward: 



T{1) 



\t{a)P 



Al + B , 



(3.7) 



where A and B are integration constants that will be determined by the physical bound- 
ary conditions associated with the temperature. We may identify in the steady states 
three different microscopic/macroscopic length (or time) scales. They result from the 
dimensionless constants 7, T{a) and A. Hence, we call them (reduced) 'reference scales'. 

In effect, ^/ri*{a) = du/dl is of the order of Vhjvy. (or equivalently, of the order of 
XrIL) and f(a) = d'^Tldt^ is of the order of {vhlvrf (i.e., of the order of {Kl^f)- 
Therefore, from (3.3 3.4|, |7|/?7*(a) and |r(a)|^/^ set two constant scales for the two 
(reduced) spatial gradients in the problem. As the reference point for computing i/^, 
we pick the one where the temperature has the lowest value, since for this choice a 
small djdl implies small djdy (we recall that djdy = yjTr/Td/dl ) . The third reference 
scale is determine d by A, which yields the reduced temperature gradient imposed at the 
boundaries: from (3.7 1 it is straightforward that A = {T{+L/2) — T{—L/2))/L, where 
±I//2 = l{y = ±h/2); therefore, A is of the order of Vh/vr- 

Thus, our three reference scales may be defined in the following way: t^ = I7I (the scale 
defined by shearing m), Tcs = |r(a)|^/^ (the scale define by collisional cooling and viscous 
heating balance) and tt = A (the scale defined by the wall temperature difference). Their 
nature and the way they affect the behaviour of the system is qualitatively different 
in each case. Two of them {ts, tt) are free; i.e., they are controlled by the boundary 
conditions. The remaining reference scale {tcs) is determined jointly by the nature of 
the granular gas (its degree of inelasticity, and thus this scale is not free like the other 
two) and the viscous heating. Once the boundary conditions are given, the three scales 
take, for the value of a of the granular gas, numerical values. We will pick the maximum 
among these three values as our definition of Knudsen number, denoted by Kn. Due to 
the existence of the inelastic cooling term, Kn cannot be made arbitrarily low as in gases 
with elastic collisions. Thus, if we look for a given type of fiow in a granular gas, there is 
a minimum value Kn„i > below which this type flow is not possible. We analyze this 
in more detail in S|4] 

f rj*{a) has little influence since 77* (a) is close t o 1 for the wh ole range of values of the 
coefficient of restitution we consider in this work (see Brey et aZ.||1998j). 
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3.3. Properties of steady temperature profiles (with or without shear) 

The change of variable we use does not change the sign of the first derivative: dT/dl has 
the same sign as dT/dy, and also dT/dl = if and only if dT/dy = 0. We must analyze 
more carefully the second derivatives: 



92f 



T 



-1/2. 



dl 



T 



-1/2 



dT 
~dl 



= T- 



, d^T 



T- 




with 



= r"2$(a) 



(3.8) 



*(a) = Bt{a) - -A^ 



(3.9) 



Relations (3.8 3.9 1 imply that, for finite temperature, the curvature of T{y) is either 
non-zero V Z,y, if <I>(a) ^ 0, or zero V Z, y if <&(«) = 0. Furthermore, the sign of <I>(a) 
yields the sign of the curvature of T(y). Thus, we call $(a), defined in (3.9 1, the temper- 



ature curvature parameter. As we see, the temperature profile curvature is determined 
by the cooling and viscous heating rate balance, through T{a), and also the boundary 
conditions for the temperature, through the constants A, B. Thus, because of the tem- 
perature dependence of the hydrodynamic coefficients, a net thermal flux can reduce or 
even reverse the effects of the heating and cooling balance on the temperature profile. 

Since d/dl = vTd/dy, we can write 



r5[ 

dy 



f{l) dl = y-yo 



(3.10) 



where yo is an integration constant, to be determined by boundary conditions. 



Integration of (3.101 involves the square root of the temperature, and for this reason, 



simpler expressions result if we first rewrite T{1) as a function of terms of the form 
(l — Iq)'^. The equality (3.7l may be rewritten, for T{a) ^ 0, as 



f([) = Lo \±ll(l-h)'^ 



(3.11) 



Given that the granular temperature cannot be negative, we have the following possible 
comb inatio ns, depending on the sign of Lq: for Lg > both upper and lower signs before 
l\ in (3.11 1 are possible, with the lower sign restricted to the interval l\{l — lo)"^ < 1; and 



for Lq < 0, only the lower (negative) sign is possible, and the solution is restricted to the 
interval ll{l — Iq)"^ ^ 1- The sets of constants {T{a),A, B} and {Lq, Ic, lo} are related as 
follows 



Ln = B- 



A^ 



$(a) 
2f (a) f (a) 



±1.^ 



T{af_ 
2$(a) ' 



k = - 



T{a) 



(3.12) 



Notice the special case ^{a) — {A^ — 2T{a)B), for which the relation (3.12 1 between 
the sets of constants {T{a),A,B} and {Lo,lc,lo} is indeterminate. In fact, the case 
A'^ = 2T{a)B corresponds necessarily to the less general form 



T{1) ^ Loiil - loi)' 



(3.13) 
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which is fulfilled if and only if A^ = 2T{a)B (i.e., <i>(a) = 0), and the relation between 
the constants is 



Lni — 



T{a) 



Icii — lo — —- 



A 



T{a) 



(3.14) 



Thus, a linear T{y) profile (because <i>(a) — 0) necessarily has the form (3.13) in the 
variable /. The existence of linear temperature profiles in granular gases heated from the 
boundaries has not been reported previously, except for a special case in a non-sheared 
granular gas (Brey et aL|2001 1. Note als o that , since the granular temperature cannot be 
negative, the lin ear temperature profile (3.131 is only possible for Lq; > 0; i.e., T{a) > 0. 



this back in (3.131 we obtain 



Inserting (3.131 in the integral in (3.101 we obtain y ^ yo ± LJ 1(1 — 21qi)/2, and using 



1/2? 



f{y) ^ ±2Ll{\y - yo) + Loilli , (3.15) 

where the plus and minus signs stand for I — Iq > and I ~ Iq < respectively. If we use 



this relation, it is straightforward, with the aid of (3.141 and (3.6 1, to find that the slope 
of the temperature profile fulfills 



dT 



dTTr 



pcrf'^^i) dy pa'^'^~^^ dy A^ 

32(d-l)7r(''-i)(l-a2) 



± 



(d-l)7r('^-i) 



(3.16) 

, 1/2 



d(d+2)2r(d/2)2(K*(a) -^*(a)) V'^ + 2/ dT{d/2)^{K*{a) - fi* (a)) ^ I 
Notice that this relation includes, as a particular case (7 = 0), the linear temperature 



described previously by Brey et al. (20011. 



In the case of a sheared system (7^0), once ( 3.10| ) is integrated the solution of the 
velocity profile is obtained using (3.3) 



77* (a) 



l + C. 



(3.17) 



where C is a constant of integration, that in the following we set equal to zero by a 
Galilean transformation. 

3.4. General classification of steady hydrodynamic profiles (with or without shear) 

We have found the following steady profiles, depending on the values of r(a), <i>(a) and 
the constants A and B: 

• Viscous heating greater than coUisional cooling {T{a) < 0, only with shear): 
o Temperature profile T{y) with negative curvature {^(a) < 0): 



T(0-Lo l-lUl-lo)' 



(3.18) 



I ri/2 1 



sm-\k{l-lo)) + lc{i-la){'^ - c^i- l^ff' 
• Viscous heating equal to coUisional cooling (r(Q;) — 0, only with shear). 
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o Temperature profile T(y) with negative curvature, but linear T{ux) (*&(q;) < 
and A 7^0): 



f{i) = Ai+B, y = yo 



'SA 



[b + Apj 



3/2 



(3.19) 



o Uniform shear flow (USF): constant temperature, linear velocity profile ($(a) = 
and A = {), B^Q): 

f{i) = B, y^v^ + B^'H. (3.20) 

• CoUisional cooling greater than viscous heating {p{a) > 0, with or without shear), 
o Temperature profile T(y) with negative curvature (<&(«) < 0), but positive curva- 
ture for r(/), T{u^): 



f{l)^L„ l-ll{l-laf 



(3.21) 



y = yo + 



\Lo\'^' 



{I - lo)^m - kY - 1 - log \21S - lo) + 2^ll{l - loY - 11 /I, 
o Linear temperature profile T{y) (<&(a) = 0): 

f{i)^Lai{i-lo)\ (3.22) 



^1/2 
y^y^±^l{l^2l^). 

o Temperature profile T(j/) with positive curvature (<&(q;) > 0): 

f(0-Lo[l + /2(f-?o)2l 



(3.23) 



y^yo + 



L, 



1/2 



1 

2L 



sinh-\l4- lo)) + Ui-lo)il + c^{l- lo)^^^ 



Flo w veloc itv profiles are obtained by making the substitution I = (ri*{a)/j)ux in 
(3.18l-(3.23l. In the case of the new linear temperature profile for T{a) > 0, among 



the two analytically possible (symmetric) solutions of the velocity profile resulting from 
(3.231, using (3.3 1 we find only one is physically possible, the one fulfilling sg{dux/dy) = 
sg{'y/r]*{a)). Note that the case r(a),$(a) = corresponds to the uniform shear flow 
(USF), which has been extensively studied ( [Campbell |1989[ ), but notice also that the 
USF that is not the only possible case with T{a) = (i.e., where viscous heati ng and 
coUisional cooling are balanced). There is also the case T{a) = 0, $(a) < 0. From (3.191, 
the range of flow velocity for this new case is limited to Ar]*{a)ux/^+B > 0, as T{y) > 0. 
The solutions we describe can be applied to speciflc boundary conditions in a straight- 
forward way. As an example, consider the boundary momentum and energy balance con- 



ditions analyzed by Alam & Nott ( 1998 1 The boundary conditions should be expressed 



Alam fc Nott|l998 1 



in the scaled variable I. The momentum balance conditions (eq. 1 in 

would flx the shear rate 7 and the energy balance conditions (eq. 2 in [Alam fc Nott 

19981, together with equation (3.7 1, would fix the constants A,B (and thus, also <I>(a 



if equation (3.9 1 is taken into account). This would determine, in the low density limit, 
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Figure 2. The values Kn_ (solid line), Kno (dotted line) and Kn+ (dashed line): (a) for a 2D 
granular gas (disks) and (6) for a 3D granular gas. 



the type of flow that would occur for given values of those boundary conditions. Either 
for these or other types of boundary conditions, the resulting type of profile should be 
described by one of the equations in (3.18)-(3.24); i.e., our classification applies for all 



types of boundary conditions, and is equivalent, when restricted to symmetric boundary 
conditions and for the low density limit, to the description of the base steady states by 



Alam & Nott (19981 and others 



It is useful to compare the sequence of steady states (3.18)-(3.24l with the picture that 
emerges from prior work. For instance, in Jenkins fc Savage j l_983[). the change of sign of 
r(a) (or, analogously. A, in the notation of Jenkins fc Savage||1983 ') implies a change of 
sign in the curvature of T(y) (the cases of negative and positive curvature in the bulk are 
also referred respectively as 'sink' and 'source' walls by |Alam fc Nott|[T998l see figures 
1 and 2 in their work). However, we see that the temperature curvature (i.e., sign of 
$(a)) and the sign of T(a) do not always coincide for asymmetric boundary conditions. 
We would expect this to have consequences on the still unexplored stability criteria of 
the new steady base states we have found. Of special interest may be the asymmetric 
counterpart of the USF (since the USF has very extensively s tudied): i.e., the new case 



with r(a) = 0, (i.e.; 'adiabatic' walls in the nomenclature of Alam & Nott 1998), but 
$(a) < 0, as described by equations (3.19). 



4. Comparison ^vith simulations and discussion 

as defined in 
Kn 



for fiows with r(a) ^ and T{a) > 0, respectively. It is also of interest to analyze 



We can calculate Kn^, as defined in § |3.2| from ( |3.5[|3.6| . Let Kn_ and Kn_|- stand for 

lalyz 
Kum for the case r(a) > with no shear (7 = 0), which we denote as Kng. In Fig. p^ we 
show Knm{a) up to the transition region from Navier-Stokes to Burnett order, usually 
taken in the range of Kn = 0.1 ( Agarwal et a?.|[2001 l. These plots provide values of the 
minimum degree of elasticity needed to allow a Navier-Stokes regime, depending on the 
type of steady base state. From figure l2] it is evident that Kn_, and thus the Knudsen 
number for USF (Kn_ indicates the value of Kn for which T{a) = 0), is always the 
largest of the three and the most rapidly growing. Other representative values, outside 
the range in the graph, are: Kn_ — 0.337 for a = 0.9, and Kn_ == 0.527 for a — 0.7. 
Nevertheless, for a > 0.995 the system has Kn_ ^ 1, and thus, there are uniform shear 
flows that may be described by the Navier-Stokes hydrodynamics. Another interesting 
result in Fig.[2]is that the sheared states are less restrictive compared to the non-sheared 
ones (as Kng > Kn_|_, systematically). Also, we see that Navier-Stokes hydrodynamics 
does not guarantee a correct description of steady shear flows in the granular gas for 
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Figure 3. The reduced normal stress differences 1 — Pxx/Pyy{^), f — Pzz/Pyy{/^) for uniform 
shear flow (r(a) = 0,A = 0,B ^ 0) in a. 3D granular gas are almost zero for a > 0.995, 
approaching the Navier-Stokes regime. DSMC simulations are performed in the dilute limit. 



a < 0.95 (although this does not mean that NS description of the base steady states will 
completely fail for a < 0.95). 

In order to qualitatively asses the reliability of the interpretation of results in Fig. [2] we 
have performed computer simulations of the USF based on the direct simulation Monte 
Carlo (DSMC) method (as implemented by Montanero et a?. | |1999 ). Concretely, we look 
at the normal stress differences, a non-linear effect (beyond NS) that appears in the 
problem for high enough Knudsen number ( Sela fc Goldhirsch||1998 1 . We plot in figure IS] 
the values of 1 — Pu/Pyy (with i = x,z) for the USF in a low density granular gas. The 
values of 1 — Pa/Pyy are very close to zero (within 1 %) for a > 0.995 but 1 — Pxx/Pyy 
starts increasing significantly (i.e. at least a departure from NS in one aspect is observed) 
for a < 0.995 (see [Montanero et al. (19991 for more simulation data in a smooth hard 
sphere granular gas and Campbell ( 1989 1 for a rough sphere gas) . The results in figure IS] 
are in agreement with the range for which Kn_ (figure p| becomes non-negligible. 

A more severe test is to compare directly the hydrodynamic steady profiles resulting 
from DSMC and MD simulations with the theoretical predictions of our study, particu- 
larly for non-homogeneous states. We have carried out DSMC simulations of hard spheres 



as implemented by Vega Reyes et al. { 2008 1 in their DSMC results for the Boltzmann 



equation. The MD simulations in this work use a standard event-driven algorithm with 
the hard sphere coUisional model (more details on methodology can be found elsewhere, 
for instance in the book by Rapaport|2004 l. In both MD and DSMC simulations we use 
'stochastic' boundary conditions ( Rapaport|2004 1 ; i.e., the heating walls act by changing 
the colliding particle velocity components so that new random values are assigned to 
them according to gaussian (for the components parallel to the wall) and Rayleigh (for 
the normal component) distribution functions with the same temperature as the pertain- 
ing wall (this type boundary conditions are commonly used in the field, see for instance 



Calvin et al. 2007). If there is shear, an a;-component constant value, characteristic of 
each wall, is added to the velocity of the colliding particle ( Vega Reyes et aZ.|2008 l. Thus, 
the input for the boundary conditions in the simulations are the fluid temperature and 
flow velocity next to the walls, and 7, T{a), and <i>(a) are extracted from steady state 
simulation data, with the aid of relations defined in §3|f[ The simulations are performed 
in dilute systems (a packing fraction 1/ ~ 7 x 10~^ is used in MD simulations) and the 
number of particles is of the order of iV == 2 x 10^ for DSMC and N = 40000 for MD. 

We discuss first the new steady states with linear temperature profile: linear T{u) 
(i.e., with f (a) = and 4>(a) < 0) and linear f{y) (i.e., with f (a) > and ^{a) = 



f If the simulation/experimental 

e.p = {d^f/dul)i^/r,'-)l,p, with (7/»7*)e. 



T{y) and 
= f^/^dujdy. 



Hy) 



measured then: 
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Figure 4. Comparison between theoretical (line), DSMC (a) and MD simulations (x) temper- 
ature profiles for two new steady states, (a) A state of the type r(Q:) = like USF but with 
$ < 0, for which theoretically T'(/), T(u) are non-uniform and linear. (6) a state with T(a) > 
and "1> = and for which theoretically T{y) is linear. In both cases a = 0.99, d = 3. In (a) the 
Navier-Stokes theory predicts 7 = 0.117, which compares (for A = 0.137) well with simulations, 
from which we measured 7 = 0.114, A = 0.149 in DSMC (2.6% and 8 .6% d iffere nce re spec- 
tively; A values are here extracted from the slope of T(m), see equations (3.191 and (3.17 1 ) and 
7 = 0.122, A = 0.140 in MD (4.6% and 4.1% difference respectively). InJT), with 7^0.091, 
the theoretical prediction for T{a) is F = 1.101 x 10~^, which compares well with the DSMC 
value f = 1.12 x 10"^ (1.7% difference) and the MD value f = 1.03 x 10"^ (6.4% difference). 



0). Figure I4] shows the results for a = 0.99 (a sufficiently high to allow Navier-Stokes 
hydrodynamics, see figure [2]). The agreement between the theoretical and the DSMC 
profiles is excellent for both cases. The deviation from a linear profile for DSMC is 
just 0.54% and 0.73% respectively (and the MD data are similarly close to linearity as 
well, see figure H caption for more details). We have performed a series of simulations 
with decr easing s hear rate and have found all the profiles corresponding to the theoretical 
solutions ( 3.18| -(3.23) (see supplementary material) . For r(a) < 0, a state with <i>(a) < 



is found. Decreasing the shear rate, we reach a flow with r(Q;) — 0, <i>(a) < (or the 
USF, with r(Q;) = 0, <i>(a) = if the fluid near both walls is at the same temperature). If 
we further decrease shear rate we find states with T{a) > with, consecutively, ^{a) < 
0,<I>(a) — 0, <i>(Q;) > (or only ^{a) > if the fluid near both walls is at the same 
temperature If] Furthermore, we found the same sequence of steady states for simulations 
far from the Navier-Stokes order, which is an indication that the qualitative behavior 
described here probably persists beyond the range of small Kn. This is perhaps not 



surprising since Tij et al. (20011 have shown that for a granular gas (and also, for a 



granular impurity, see Vega Reyes et al. 2008), the Couette flows obey hydrodynamic 



equations that share the same structure with those at Navier-Stokes order, albeit with 
new nonlinear transport coefficients and r(a), dependent on a constant shear rate, that 
would have to be calculated ( |Tij et a?.||2001| [Santos e^"ar]|2004[ [Vega Reyes et aZ.||2008 1 . 
MD simulations results in figure l5l show another type of profile not reported (to our 
knowledge) in the literature on granular gases in contact with thermal walls. In this case 
there is no shear in the system (7 = 0) and the temperature profile shows the nega- 
tive curvature characteristic of an elastic gas subjected to two thermal walls at different 
temperatures. As in the linear temperature profiles, there is good agreement between 
simulation and theory. The boundary conditions we used for the MD simulations are 



very common in the literature (Calvin et al. 20071, but one could use in the MD code 



\ We refer the reader to the supplementar y mat er ial ac companying this article for figures of 
all different types of flows, as extracted from (3.18l-(3.23l. 
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Figure 5. Temperature profile with negative curvature (elastic gas- like) for a granular gas 
(a = 0.999) heated from the walls and no shear (7 = 0). This type of profile has been not reported 
before in granular gases. The line corresponds to the theoretical profile with $ = —0.0321 and 
points correspond to MD simulation data, with $ — —0.0318 (1% difference; the theoretical 
profile was adjusted to make t he te mperature boundary values coincide, and $ in the simulation 
is obtained using the relation (3.8l). 



boundary conditions closer to experimental boundary conditions; for instance, bumpy 
moving walls to produce shear ( Lun|199"6 l, or vertically vibrating walls to heat the gran- 
ular gas ( Vega Reyes fc Urbach|2008 1. It is reasonable to expect good agreement between 
our theory and the simulation results with these more realistic boundary conditions, since 



it has been found that this agreement exists in experiments (|Yang, Huan, Candela, Mair 
fc Walsworth||2002 1 for the steady states derived previously. 



The physical origin of the sometimes surprising behavior of the temperature profile is 
the temperature dependence of the hydrodynamic transport coefficients, combined with 
the effects of constant pressure. Consider, for example, the case of a linear temperature 
profile when collisional cooling exceeds viscous heating. If the transport coefficients were 
approximately constant, the heat ffux would be constant across the system, and there 
would be no source to replace the energy lost to collisional cooling. However, in the 
hydrodynamic equations the heat ffux is actually proportional to vTdT/dy, and therefore 
is not in general constant for a linear temperature profile. In the case of the linear profile, 
the decrease in heat flux from the hot side of the system to the cold side compensates for 
the collisional cooling. A perhaps surprising result is that if this balance is satisfied at 
one point in the fiow, it is satisfied everywhere (or, more generally, that the curvature of 
the temperature profile is constant, even though the temperature and density are not). 
This is because the hydrodynamic fields scale in the same way as the pressure, which is 
constant across the system. 

The results for DSMC and MD simulations supply information at two different levels. 
DSMC results show that our theory, resulting from an expansion of the solution of the 
kinetic equation, approaches in effect the exact solution of the kinetic equation (yielded 
by DSMC). MD simulations help to check that the description resulting from the kinetic 
equation approaches a 'real' system. The states that we have described in molecular 
dynamics simulations should be observable, in principle, in experiments. The temperature 
and/or heat ffux of the sheared granular gas at the boundary is determined by the 
physical properties of the shearing walls (roughness, friction, vibration, etc.). Starting, for 
example, from symmetric relatively smooth non-vibrating walls, and then systematically 
increasing the roughness and/or vibration amplitude of one wall, the ffow states we have 
identiffed can be explored. 

The unidimensional steady states analyzed here are the starting point for more complex 
studies on granular ffuid ffow (|Batchelor||1967|). For example, the fact that the linear 
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T{u) shows exactly the same transport properties and rheology as the USF in the regime 



beyond Navier-Stokes order (Vega Reyes et al. 2008 Santos et al. 2009), will allow the 



measurement of heat flux transport coefficients directly from a numerical solution of 
the Boltzmann equation in steady state, something that was not possible previously 
(|Garz6 2006). Preliminary DSMC results confirm that this state also appears beyond 
NS order IfJ and that the (non-linear) shear viscosity values (from DSMC) for both the 
USF and the new states show nearly identical values (Santos, Vega Reyes &; Garzo 



[2009 j ) . Another possible application on which we are working is the Brazil-nut segregation 
problem, since the input of the solutions of the profiles is needed in order to determine 
the segregation criteria (Garzo fc Vega Reyes|2009 1. Additionally, results in figure Is] show 



that the positive curvature caused by coUisional cooling can be reversed if the thermal 
walls are at sufficiently different temperatures. Thus, the results allow the extension 



of the study of deviations from NS hydrodynamics for a heated granular gas (Hrenya, 



Galvin fc Wildman_2008) (no ffow velocity) to the three classes of steady states that will 



appear, depending on the wall temperature difference, as described by (3.21|-(3.24l. It 
would also be interesting to study the effect that this heat ffow gradient sign reversal 
mechanism would have on hydrodynamic stability problems applied to granular gases, 
such as the buoyancy instability mechanism in a ffuid layer between two walls at different 
temperatures (Rayleigh-Benard convection). 
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